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NON-LINEAR BEHAVIOR OF 
FIBER COMPOSITE LAMINATES 

by Zvi Hashin, Debal Bagchi and B. Walter Rosen 
Materials Sciences Corporation 

SUMMARY 

The non-linear stress-strain behavior of fiber composite 
laminates has been analyzed to define the relationship between 
laminate behavior and the non-linear stress-strain characteris- 
tics of unidirectional composites. The resulting analysis has 
been programmed to yield an efficient computerized design and 
analysis tool. 

The approach utilized herein was to adopt a Ramberg- 
Osgood representation of the non-linear stress-strain behavior 
and to utilize deformation theory as an adequate representation 
of the material nonlinearities. The problem was viewed on two 
levels. First, the relationship between the constituent proper- 
ties and the stress-strain response of a unidirectional fiber 
composite material was studied. For this problem, the primary 
attention was directed toward axial shear behavior, and an ex- 
pression was established relating the composite average-stress/ 
average-strain curve to the fiber moduli and the matrix non- 
linear stress-strain curve. Second level of approach is to treat 
the interelationship between the properties of the unidirectional 
layers and those of the laminate. For this case, the starting 
point is a non-linear stress-strain curve for transverse 
stress and for axial shear and a linear stress-strain relation 
for stress in the fiber direction. The non-linear lamina 
stress-strain curves can be modeled by proper selection of the 
Ramberg-Osgood parameters. In the present study, with this as a 
starting point, an interaction expression was formulated to ac- 
count for simultaneous application of axial shear and transverse 
stress. 

A laminate having an arbitrary number of oriented layers 
and subjected to a general state of membrane stress was treated. 
Parametric results and comparison with experimental data and 
prior theoretical results are presented. 
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INTRODUCTION 


A basic requirement for the engineer designing with fiber 
composite materials is a definition of the stiffness and strength 
of these materials under a variety of loading conditions, includ- 
ing cases for which experimental materials properties data are 
not available. For this purpose, it is necessary that he have at 
his disposal reasonably accurate procedures to predict these 
mechanical properties. Existing analyses can predict the elas- 
tic behavior of a laminated composite quite well when the elas- 
tic properties of the unidirectional materials from which it 
is made are known. However, the situation has been much more 
complicated and much less satisfactory with regard to the in- 
elastic stiffness and strength of a laminate. The present pro- 
gram was undertaken to develop a computerized analysis of the 
inelastic behavior of fiber composite laminates which could be 
used as a design tool. The results of this study and compari- 
sons of these results with experimental data are presented in 
this report. 

It is essential to recognize that the utilization of fiber 
composite materails in structural design involves the incor- 
poration of material design into the structural design process. 
This is illustrated clearly by the fact that the gross mate- 
rial properties of a fiber composite laminate change when any 
change is made in the laminate ply orientations. Even when 
the designer considers a material formed from a particular com- 
bination of fiber and matrix materials, there remains a large 
number of geometric variables associated with the laminate de- 
sign. Thus, in the preliminary design phase, experimental mate- 
rial properties data will generally be too limited. In the 
case of elastic properties, sufficient capability to synthesize 
the necessary properties exists. This procedure generally starts 
with the definition of the elastic properties of unidirectional 
fiber composite materials. These can, of course, be determined 
experimentally. Also, when such data are not available, they 
can be estimated using a variety of analytical techniques. These 
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latter are generally referred to as micromechanics analyses. 

For example, a set of relatively simple relations for predict- 
ing the moduli of unidirectional reinforced composites are 
presented in [1], Alternate micromechanics approaches are 
described in [2] to [4] . A review of these methods .is pre- 
sented in [5]. With these properties available, it is assumed 
that the individual laminae are homogeneous and anisotropic. 

A laminate analysis is carried out in a straight forward 
fashion following methods originally developed for such mate- 
rials as plywood, and more recently extended to the more 
general cases associated with fiber composite laminates (e.g., 
[6] to [8]). 

However, contemporary fiber composite materials generally 
consist of elastic brittle fibers such as glass , boron or 
graphite in relatively soft matrix materials such as epoxy or 
aluminum. For these matrix materials it is reasonable to an- 
ticipate that at a certain loading state the matrix will begin 
to exhibit inelastic effects. This results in non-linear re- 
lations between structural loads and deformations. These in- 
elastic effects can, of course, be expected to have a signifi- 
cant effect upon failure of the laminate. It is quite clear 
that adequate definition of these failure conditions are 
essential to achieve structural designs of high reliability. 

In the present study, a non-linear laminate analysis has 
been developed which can provide realistic assessments of the 
stresses and strains in the various laminae and of the inelas- 
tic stiffnesses of the laminate at any stress level. This in- 
formation can be used for assessment of such effects as struc- 
tural stability or structural stress distributions. The stress 
distributions in the laminae and the laminates can also be 
utilized for the development of more realistic failure criteria 

Inelastic matrix behavior can be classified broadly as 
either time dependent or time independent. Time dependent be- 
havior is called viscoelastic if linear and creep if non-linear 
Polymeric matrices such as epoxy do exhibit such behavior. In 
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the case of metallic matrix materials, such as aluminum, time 
dependent effects are generally negligible unless elevated 
temperature conditions are considered. The present study is 
concerned with time independent non-linear matrix behavior 
which is of significance for both polymeric and metallic matri- 
ces. Throughout this paper the expression "inelastic" is used 
to describe this time independent mechanical behavior. The 
method of approach to these problems is similar to that of the 
elastic analysis. Thus, it is necessary to determine, first, 
the inelastic properties of the- unidirectional fiber composite 
materials. This can be done experimentally or by micromechan- 
ics methods. Given this information, a method to determine 
stresses and strains in an inelastic laminate is then devised. 
The problem is complicated by the fact that the inelastic 
stress-strain relations are non-linear. 

A limited number of pertinent investigations can be found 
in the literature. Hill [4] considered, in approximate 
fashion, a limited aspect of inelastic behavior of a uniaxially 
reinforced material: the case of stress in fiber direction 
combined with isotropic transverse stress. Petit and Waddoups 
[6] devised an incremental method for laminate analysis in 
which it was assumed that in single laminae there is no inter- 
action of stress components in different directions as far as 
lamina deformation is concerned. This assumption is restric- 
tive, and also their incremental laminate analysis scheme is 
unduly complicated. Adams [7] used a finite element technique 
for numerical analysis of unidirectional materials in the form 
of periodic fiber arrays under conditions of plane strain. 

Huang [8] gave an approximate analysis for transverse inelas- 
tic behavior for a unidirectional material in plane strain, 
but it is diffucult to assess the validity of the approxima- 
tions introduced. 

A detailed analysis of the inelastic laminate problem 
has been given by Foye and Baker [9] . Using finite element 
methods, they computed the inelastic effective properties of 
unidirectional rectangular and square arrays of elastic fibers 
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in inelastic matrix. These properties were then used in an 
inelastic laminate analysis. The analysis is based on 
incremental plasticity theory and is, unfortunately, very 
complicated and requires a great deal of computer time. The 
results obtained are, however, of great importance for com- 
parison with results predicted by more simplified theories, 
such as the one which will be given in the present work. 

The body of this report is divided into four major sec- 
tions. In the first, consideration is given to the behavior 
of unidirectional fiber composite materials. This requires: 
a definition of the appropriate form of the inelastic stress- 
strain relations; some consideration of the relationship be- 
tween composite properties and constituent properties; and a 
definition of the appropriate form of the interaction between 
various stress components. The basic objective in this phase 
of the report was to define appropriate constitutive relations 
for the individual lamina which can be used in the non-linear 
laminate analysis. Further, there is a desire to gain some 
insight into the influence of the particular constituent 
properties upon the lamina stress-strain relations. In this 
phase of the study, it is found useful to characterize the 
unidirectional material with the aid of Ramberg-Osgood stress- 
strain relations. 

In the next section of the report, the analysis of the 
inelastic behavior of laminates is described. Here, a pro- 
cedure for incorporating the non-linear constitutive relations 
into an analysis which defines the state of stress in the in- 
dividual laminae under an arbitrary set of external loads, 
is defined. Analyses are developed for the case of symetric 
laminates subjected to membrane loading. The equations which 
are developed uniquely define the desired laminate internal 
average stress distribution under a given set of membrane loads. 
Governing equations, however, are non-linear and require numeri- 
cal solution procedures. An efficient algorithm has been de- 
fined which enables computer solution to be achieved for arbitrary 
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laminates at minimal cost. The solution is obtained by appli- 
cation of the Newton-Raphson method. 

In the final section, the computerized analysis which has 
been developed is applied to series of problems. The first 
group presents comparisons with ‘various analytical results from 
the more complex analyses of Ref. [6] and [9]. The second group 
of numerical results presents comparisons between theoretical 
results from the present model and available experimental data. 
The third group of results provides several parametric studies 
to gain insight into those factors which contribute signifi- 
cantly to the non-linear behavior of fiber composite laminates. 
Also, computations have been made to provide a preliminary 
assessment of combined load effects including comparisons with 
limited experimental data. 

Details of the various analytical developements , as well 
as descriptions of the computer program, are presented in 
appendices to the report. 

The principal result of the present program is a computer 
program which provides a simple engineering tool which can be 
used for the parametric study of the influence of material prop- 
erties upon laminate performance. This laminate analysis capa- 
bility can be used by the structural designer to define design 
allowable stresses and to aid in the selection of fiber com- 
posite materials for structural applications. A comparison of 
the present results with the limited amount of available experi- 
mental data shows good agreement. There are, however certain 
cases in which the agreement is not good, par ticularly as the 
laminate loading approaches failure. The results of the present 
analytical method agree well with the results for those prob- 
lems for which more exact and more complex analytical results 
exist . 
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2 . 


NON-LINEAR STRESS-STRAIN RELATIONS OF UNIAXIAL 
FIBER REINFORCED MATERIALS 


2 . 1 General Form of Stress-Strain Relations 

An effective stress-strain relation of a composite mate- 
rial is defined as a relation between average stress and 

average strain i^j. Here and in the following latin indices 
range over 1, 2, and 3. If the composite is elastic the 
general effective stress-strain relation takes the form 


a ij - C *ijkl 


'kl 


( 2 . 1 . 1 ) 


where are the effective elastic moduli which are material 

constants and are thus independent of stress or strain. Thus, 
( 2 . 1 . 1 ) is a linear relation between average stress and strain. 

If the composite is subject to symmetries the form of 
(2.1.1) simplifies. For a uniaxial FRM the most important 
cases of symmetry are transverse isotropy, around fiber direc- 
tion, and square array (square symmetry). In these cases the 
stress-strain relations ( 2 . 1 , 1 ) for transverse isotropy assume 
the form: 


a l l 

II 

n 

* 

£.11 

+ c* 12 

£ 2 2 

+ 

C * 1 2 

£ 3 3 

0 2 2 

= c\ 2 

e i i 

+ c* 

^22 

E 2 2 

+ 

c*. 

^ 2 3 

F 3 3 

O 3 3 

II 

O 

* 

IO 

£ l l 

+ C*2 3 

£ 2 2 

+ 

C*2 2 

£ 3 3 


o 1 2 — 2C* 4 4 e 1 2 


02 3 — 2C* 5 5 £23 

Oil = 2C* 4 4 £ 3 1 

and 

C*5S= (C* 2 2~C* 2 3 ) /2 


(2.1.3) 


In (2. 1.2-3) 1 indicates direction and 2, 3 perpendicular direc- 
tions transverse to 1 . 

In the event of inelastic matrix and elastic fibers, the 
situation is much more complicated since the stress-strain 



relation are nonlinearity and history dependent. In no case is 
stress proportional to strain so that superposition of effects 
is not valid, and in order to determine current strain it is 
not sufficient to know current stress but it is necessary to 
know precisely the variation of stress which preceded its cur- 
rent value. Thus, for a material in a known state of combined 
shear and uniaxial tension, the state of strain is different 
if: (a) tension is first applied and then the shear, (b) shear 

is first applied and then the tension- (c) tension and shear 
are applied simultaneously. For this reason stress-strain re- 
lations must be presented in incremental form. That is, strain 
increment is related to stress and stress increment. This com- 
plicates matters enormously. However, it is known that in 
the case of proportional loading, that is, all stresses at a 
point grow simultaneously in a fixed ratio to one another, 
incremental theory can be integrated into the much simpler 
total or deformation theory for which current strain is com- 
pletely determined by current stress. 

Deformation theories have a wider range of validity than 
proportional loading. Comparison of numerous detailed solutions 
carried out both incrementally and by much simpler deformation 
theory show surprising agreement in many cases, and Budiansky 
[10] has shown that deformation theory can also be valid for 
"neighboring" loading paths. 

In the present work, we are concerned with composites 
which are subjected to some external load. If it is supposed 
that the various external load components grow proportionally, 
this does not necessarily imply that the components of stress 
at a typical internal point also grow proportionally. It is, 
however, felt that the manner of growth of these internal 
stress components cannot deviate severely from proportional 
loading if external loading is proportional. Consequently, 
deformation type stress-strain relations are assumed for the 
matrix. 

This assumption results in considerable simplification. 

It will be seen that it yields results which are extremely 
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close to the ones obtained in [ 9 ] on the basis of the much 
more complicated incremental theory. 

It is shown in Appendix A that for elastic fibers and 
an inelastic matrix described by deformation type theory, the 
effective stress-strain relations for a transversely isotropic 
or square symmetric FRM are: 


£ 1 1 

= S 1 1 

On 

+ 

S 1 2 

0 2 2 

+ 

S 1 2 

o 3 3 


£22 

= S 1 2 

Oil 

+ 

S 2 2 

0 2 2 

+ 

S2 3 

0 3 3 


£33 

- s 1 2 

On 

+ 

S2 3 

0 2 2 

+ 

S 2 2 

0 3 3 

( 2 . 1 . 4 ) 


£ ! 2 = 2 S 4 4 3 ]2 

£23 = 2 Sss 0 2 3 
£ 1 3 = 2 s 4 4 O 1 3 

and 

Ssj= (S22-S2 b )/2 ( 2 . 1 . 5 ) 

The coefficients etc * are the effective inelas- 

tic compliances of the material and are functions of the aver- 
age stresses , or rather of certain invariants of the average 
stress tensor. 

We are here primarily concerned with thin uniaxially re- 
inforced laminae which are in a state of plane stress. Let 
x^ denote fiber direction, x 2 direction transverse to fibers 
in lamina plane, and x^ direction perpendicular to lamina. 
Figure 1 . Then the plane stress condition is expressed by: 

013 = 023 = 033=0 ( 2 . 1 . 6 ) 

Equs . ( 2 . 1 . 4 ) then assume the form: 


£ 1 1 

= S 1 1 

0 1 1 

+ S i 2 

0 2 2 

( 2 . 1 . 7 ) 

£22 

= S 1 2 

0, i 

+ S 2 2 

02 2 


£12 

= 2S 4 4 

0 1 2 
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Note that does not vanish. It is however of no interest 
for present purposes. 

The inelastic compliances in ( 2 . 1 . 7 ) are functions of 
the stresses o^, o 22 °12’ 

It is convenient to split the strains in ( 2 . 1 . 7 ) into 
elastic strains e a , and inelastic strains e e . Thus: 


e 


a$ 



+ 



( 2 . 1 . 8 ) 


where here and in the following greek indices range: over 1, 2. 
The elastic strains are recovered after unloading of the com- 
posite and are related to the stresses by elastic stress-strain 


relations . 

Thus : 




— | 

C 1 1 

= Si! 

Oil 

+ S ! 2 

a 2 2 

£22 

CM 

II 

0 1 1 

+ S22 

CM 

CM 

1 D 

_ 1 

£ 1 2 

— 2 S 4 4 

0 1 2 




where 




"A 


1 


4 G 


A 


( 2 . 1 . 10 ) 


Here is the effective Young's modulus in fiber direction, 

v,- the associated effective Poisson's ration, E m - the effec- 
A T 

tive Young's modulus transverse to fibers and G^ - axial effec- 
tive shear modulus, related to 1-2 shear. 

The inelastic, permanent, strains then have the form: 


_ 11 

11 _ 


«i. 


£ 1 1 

= S 1 1 0i 

1 + 

S 1 2 

0 2 2 

_ n 

11 __ 


11 


£22 

= S 1 2 0 

1 1+ 

S2 2 

O 2 2 

— 11 

II _ 




E 1 2 

= 2 S 4 4 0 

1 2 




where 

II II _ _ 

= 1 1 1 02 2 , 0 1 2 ) 


( 2 . 1 . 11 ) 


( 2 . 1 . 12 ) 
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In order to further simplify the stress-strain relations 
(2 . 1 . 11- . 12) , some specific features of FRM will be taken into 
account. In such materials, the fibers are by an order of mag- 
nitude stiffer than the matrix (for the case of boron and/or 
graphite fiber in an epoxy matrix the ratio of fiber to matrix 
Young's modulus can be in excess of 100). The stiffness ratio 
becomes larger in the inelastic range since the matrix loses 
stiffness (i.e., flows) while the fibers retain their stiffness. 
It is, therefore, clear that the stress in fiber direction 

is practically carried by the fibers alone, with insignificant 
matrix contribution. 

On the other hand, the transverse stress a 22 and the shear 
stress a 12 are primarily carried by the matrix with little 
fiber contribution. 

It follows that inelastic behavior of the FRM is produced 
primarily by a 22 and a^ 2 while inelastic behavior for o^ load 
can be neglected. 

The foregoing comments are summarized into two basic 
assumptions : 

_ H _ It 

(a) the inelastic strains e 22 and c^ 2 are not functions 

of 5 X1 

_ ii 

(b) the inelastic strain always vanishes. 

On the basis of these assumptions, the stress-strain re- 
lations (2.1.11-.12) simplify to: 


— II 


e i i 

= 0 



_ ll 

ii _ 



e 2 2 

= S2 2 (0 22 f 

0 1 2 ) 

02 2 

_ II 

II _ 



£ 1 2 

= 2 S 4 4 (CJ22/ 

0 1 2 ) 

0 1 2 


(2.1.13) 
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2 . 2 Plane Stress-Strain Relations in Rambera-Osaood Form 


A convenient representation of non-linear one dimensional 
stress-strain relations has been given by Ramberg and Osgood 
[11]. For uniaxial stress, for example: 




+ k (- 




( 2 . 2 . 1 ) 


where represents the elastic Young's modulus, and k, o', 
and m are three parameters to be obtained by curve fitting. The 
parameter o' is sometimes called nominal yield stress. Equa- 
tion (2.2.1) represents a family of curves with initial slope 
E^ , and monotonically decreasing slope with increasing a. The 
curves flatten out with increasing m (Fig. 2) . Without loss 
of generality (2.2.1) can be written in the form: 


e 



[1 + (^r) m_1 ] 


( 2 . 2 . 2 ) 


which will be used from now on. Similarly, a stress-strain 
curve in shear can be represented in the form: 


T 

Y “ [1 + (I_) n-1 ] (2.2.3) 

T 

y 

where is the elastic shear modulus. 

It should be emphasized that (2. 2. 2-. 3) are valid only for 
one dimensional cases. The question of the generalization to 
general states of stress and strain has no unique answer. One 
common used form is isotropic J ^ deformation theory [12] . 

Next, we consider the case of effective or macroscopic 
stress-strain relations for the special case of a uniaxially 
reinforced material in which the matrix in non-linear, with 
stress-strain relations in Ramberg-Osgood form. 

Consider, for example, the case of uniaxial average stress 
022 in direction transverse to fibers, all other average stresses 
vanish. It then follows from (2.1.7) that: 

£22 = S 22 (O 22 ) O 22 (2.2.4) 


12 . 


Similarly, if the only nonvanishing average stress is o ^ 2 
the shear stress-strain relation of the composite is: 


£ 1 2 — 2S 44 (O 1 2 ) 012 


(2.2.5) 


Evidently the inelastic effective compliances and 
are functions of the parameters of the inelastic Ramberg- 
Osgood stress-strain relations of the matrix, of the elastic 
properties of the fibers and of the internal geometry of the 
composite. Actual prediction is a very difficult problem. Such 
problems will be cnsidered in limited fashion in the next para- 
graph. 

Just as matrix stress-strain relations are represented in 
Ramberg-Osgood form, the same type of curve fitting can also 
be applied for the effective stress-strain relation of the 
composite. Thus (2. 2. 2-. 3) are written in the form: 


°22 022 u.l 

■ E 11 + > 1 

T y 


£12 = U ^ [1 + (^.) N - 1 ] 

2G A T y 


(a) 


(b) 


( 2 . 2 . 6 ) 


Where E T . is the effective transverse elastic Young's modulus 
- effective axial elastic shear modulus and o , t , M and 
N are curve fitting parameters which are in general quite 
different from the corresponding Ramberg-Osgood matrix para- 
meters . 

A question of fundamental and of practical importance is 
the form of the stress-strain relations for the case of plane 
stress, taking into account interaction among the various stress 
components. It should be noted in this repsect that (2.2.6) 
are special stress-strain relations when a 22 or ® \2 act onl Y 
by themselves. 

It is recalled that equations (2.1.13) represent the 
inelastic parts of the strains for plane stress-strain re- 
lations for FRM with stiff fibers. It is shown in Appendix B 
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that the Ramberg-Osgood form of such plane stress-strain re- 
lations is as follows : 


£n = 0 

— n _ O 2 2 

£22 “ p 


£ 12 = 


[(^) 2 

Y 

[(^) 2 

Y 


+ (^) 2 ] 

y 

+ (f^-) 2 ] 

y 


The parameters E^, , G^, a^, t , M, N in (2.2.7) are those of 
the one dimensional stress-strain relations (2.2.6) which may 
be regarded as experimentally (or perhaps theoretically) known. 

The inelastic parts of the strains are given by (2.1.9- 
.10), and the total strains are then given by adding equations 
(2.2.7) and (2.1.9) . 


Equations (2.2.7) have been compared with computed numeri- 
cal results given in [9]. Reasonable agreement was obtained. 
Comparisons for the interaction cases of transverse stress, o 22 , 
versus transverse strain, e 22 , the P resence °f axial shear 
stress, o 12 , and axial shear stress, 5 12 , versus axial shear 
strain, Y 12 / are shown in Figures 3 and 4 respectively (in both 
cases & 22^ a 12 = 8 /^) * ■‘■ s seen that the agreement is fair 

for transverse stress-strain relations (Fig. 3) and very good 
for the shear stress-strain relations (Fig. 4) . 

Figures 3 and 4 also show the stress-strain relations ob- 
tained from Eqs . (2.2.7) for one dimensional transverse tension 

o 2 2 , and axial shear, a^ 2 » respectively. 
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2.3 Axial Shear Stress-Strain Relation 


This paragraph is concerned with the problem of prediction 
of a one dimensional effective axial shear stress-strain re- 
lation of a uniaxial FRM in terms of matrix and fiber properties 
and the internal geometry of the composite. 

The main reason for concentrating on the axial shear prob- 
lem is that the inelastic effect is predominant in axial shear 
for which significant nonlinearity of the stress-strain response 
is obtained (e.g. , Figure 4) . The effect in fiber direction is 
practically non-existent as has indeed been assumed above, and 
is relatively small in transverse stress which is shown by the 
small curvature of the stress-strain relation in this case 
(e.g. , Figure 3) . 

On the basis of all this , it can indeed be assumed as 
first approximation that the nonlinearity of the uniaxial FRM 
is limited to axial shear alone. 

Consider a uniaxially reinforced lamina which is subjected 
to pure axial shear, Figure 5, on its surface. The boundary 
conditions are : 


X 3 

= + 

t/2 

a 3 1 

- O 3 2 

— O 3 3 

= 0 


x 2 

= + 

b 

a 1 2 

= 

0-2 2 

= 0 2 3 

= 0 

X 3 

= + 

a 

a i 2 

= To 

O 1 1 

= Oj 3 

= 0 


It may be shown that under such load the only nonvanishing 
average stress in the composite is : 


It would seem at first that, given the complexity of the in- 
ternal geometry of the composite, the state of stress at any 
interior matrix or fiber point is generally three dimensional. 
Surprisingly enough, however, this is not so and the only non- 
vanishing stress components in the interior of the composites 
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are the shear stresses o-^ and o^ 3 , which are moreover func- 
tions of x 2 and only. Thus, the interior state of stress 
is : 


0 1 2 

— a 1 2 

(x 2 , x 3 ) 

a j 3 

= o i 3 

(x 2 , x 3 ) 

Oil 

= a 2 2 

= a 3 3 = 


(2.3.3) 


The validity of equations (2.3.3) for the case of an elastic 
composite has been proved in [5], Their validity for the pre- 
sent much more general inelastic case will be shown elsewhere. 

The effective stress-strain relation of the composite in 
axial shear is defined by: 

T 


e 1 2 = 


O 1 2 
2G? 


2g; 


(2.3.4) 






g 

where G, is the effective secant shear modulus of the mate- 
A 

rial. The nonlinearity of the stress-strain relation is ex- 
pressed by the fact that function of the applied stress. 

. . s 

It is seen that xn order to determine G^ it xs necessary 

to compute the average shear strain for given applied 

shear stress. This is a formidable problem even with the 

simplification (2.3.3) and we shall content ourselves with a 

brief outline of its formulation. To simplify matters, the 

fibers shall be assumed to be ideally rigid relative to the 

matrix. This is a very accurate assumption for the case of 

Boron and Graphite Fibers. There is no difficulty to extend 

the formulation to the case of non-rigid elastic fibers. 

In view of (2.3.3) the problem is two dimensional and 

need only be considered in a typical x 2 , x^ section. In the 

matrix domain: 
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+ 


0 


(2.3.5 a) 
(2.3.5b) 


3o i 2 
3X2 

e i 2 = 


3o 


3x 


i 3 
3 


0 1 2 
2G 


[1 + 



n-1 


] 


- u 1 3 


2G 


[1 + (- 


n- 1 . 


T 


*0 * 
0 1 2 


+ O 


1 3 


(2.3.5c) 

(2.3.6) 


£ 1 2 
£ 1 3 


1 

2 

1 

2 


u 


1 


3 u i 
3x 2 

3u i 

3x 3 

(x 2 ' X 3 } 


(2.3.7a) 
(2. 3.7b) 
(2.3.8) 


and, u 1 = 0 at fiber/matrix interface. 

Here equ. (2.3.5) is the only surviving equilibrium 
equation, (2.3.6) are Ramberg-Osgood stress-strain relations 
for isotropic J 2 theory (2.3.7) are usual strain-displace- 
ment relations in which u 2 and u 3 do not enter since it may 
be shown that they are not functions of x^ and (2.3.8) ex- 
presses the ideal rigidity of the fibers. 

Equs . (2. 3. 5-. 8) must be solved subject to boundary 

condition (2.3.1). If this is done the strain e^ 2 i- s known 
everywhere and can be averaged to obtain G^ from (2.3.4). 

The problem is exceedingly difficult because of the non- 
linearity introduced by the stress-strain relations (2.3.6). 
There is very little hope to solve it analytically for any 
kind of fiber geometry. It should therefore be handled by 
numerical methods for fiber arrangements and fiber shapes 
of engineering interest. 

Another way to approach an analytically intractable 
problem such as the present one is by variational techniques • 
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In this fashion, approximations or bounds for quantities 
of interest are obtained by methods which are much simpler 
than bonafide solution of the problem. Such variational 
methods have been extensively used for determination of 
effective elastic moduli of FRM (e.g., [5]). 

In the course of the present work, it has been found that 
variational methods can also be used for inelastic problems 
such as the present one to obtain bounds on effective secant 
moduli. The main ingredients of the method are: 

(a) Construction of an extremum principle in terms 
of an energy integral such that the true energy 
is the minimum of the integral. 

(b) Expression of the true energy in terms of effec- 
tive secant modulus. 

(c) Establishment of admissible fields to obtain a 
value of the energy integral which is larger than 
the true energy, thus obtaining a bound for G^. 

The work involves complicated developments and deriva- 
tions which are given in Appendix C. Here only the end result 

s 

for a lower bound on will be given for a special geometry 
of FRM which is known as composite cylinder assemblage. This 
geometry has been described in detail in [1, 5] and consists 
of an assemblage of composite cylinders of variable sizes 
which are joined together so as to fill the whole volume of 
the composite. In order to fill the whole volume, composite 
cylinders vary from finite to infinitesimal size. This geom- 
etry has been used to advantage for elastic FRM to obtain 
simple expressions for effective elastic moduli which are well 
verified by experiment [1, 5]. In the present case only a 
lower bound on has been obtained for the case in which the 
exponent n in matrix stress-strain relations is n=3. 

It has been found that with this exponent and proper 
choice of t , epoxy shear stress-strain relations can be well 
described. The result for the lower bound is: 
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(2.3.9) 


g a - g a(-) 


1+c 

1-c 


1 + <i*) 2 3^3c*c 2 , tc i 
t' 3(1+0 


where 

c - volume fraction of fibers 
G - elastic (initial) matrix shear modulus 
t' - Ramberg-Osgood matrix stress parameter, and 

t q - applied shear stress. 


It follows from (2.3.4) that: 


£ 1 2 


< 


To 


2G 


s 

A(-) 


(2.3.10) 


In other words , with the lower bound on G^ an upper bound 

on e. _ variation with t is obtained. 

12 o 

If (2.3.10) is explicitly written in terms of (2.3.9) it 
assumes the form: 


£ 1 2 


2Gy-^£ 

1-c 


[1 + 



3+13c+c 2 +c 3 , 

t — 1 

3 (1+c) J 


(2.3.11) 


Recalling that for the composite cylinder assemblage 
with rigid fiber.s the axial elastic shear modulus G^ is given 
in [1, 5] as: 


G 


A 


G 


1+c 

1-c 


(2.3.12) 
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and comparing (2.3.11) with (2.2.6) with choice of exponent 
N=3 (which is the same as matrix exponent) , it is seen that: 


T 


2 

y 


i 2 

> T 

- y 


3(l+c) 3 
3+13c+c 2 +c 3 


(2.3.13) 


The prediction of (2.3.11) has been compared with numeri- 
cal results obtained in [9]. Figure 6 shows the variation of 
the right side of (2.3.11) in comparison with the results ob- 
tained in [9] for a fiber volume fraction, c=0.5. Since re- 
sults of [9] were for boron fibers in epoxy matrix, the rigid 
fiber approximation is accurately valid. It is seen that the 
results are reasonably close. It should be noted that the 
geometry of [9] is a rectangular fiber array which is quite 
different from the composite cylinder assemblage geometry. 

The results defined by (2.3.12) and (2.3.13) used in 

equation (2.2.6) yield the result plotted in non-dimensional 

form in Fig. 7. The shear strains are normalized with respect 

to the matrix elastic strain, y , at the yield stress, x : 

ye y 



It is natural to also consider the establishment of an 

g 

upper bound on G^. Unfortunately, however, this is a matter 
of formidable difficulty for the reason that inversion of 
(2.3.6) to express stresses in terms of strains leads very 
complicated expressions. Further discussion of this diffi- 
culty is given in Appendix C. 
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3. 


ANALYSIS OF NON-LINEAR LAMINATES 


3.1 Formulation 

The general problem to be investigated in the present 
chapter is as follows: given the inelastic stress-strain re- 

lations of uniaxially reinforced laminae determined theoreti- 
cally or experimentally, and a laminate composed of such 
laminae and loaded on its edges by uniformly distributed loads 
in the plate of the laminate : 

(a) What are the stresses in the various laminae? 

(b) What is the macroscopic strain response of the 
laminate to the loads? 

This problem has been extensively investigated for elastic 
laminates, and the results obtained will serve as important 
guidelines for the present much more complicated problem. It 
is therefore very helpful to first briefly review the theory 
of elastic laminates. 

Let the laminate be referred to a fixed system of coord- 
inates x^, x 2 , x^ as shown in Figure 8. This will henceforth 
be referred to as the laminate coordinate system . 

Any lamina, kth say, in the laminate will be referred to 


(k) 


(k) 


x^ where x^ 


(k) 


its material system of coordinates x. 

nrr 

is in fiber direction, perpendicular to fiber direction 

and x^ is the same as the laminate x 3 , Figure 8. The reinforce- 
ment angle 0^ is defined by: 


9 = } (x 1 , x x (k) ) = } (x 2 , x 2 (k) 


(3.1.1) 


Let it be assumed that the laminae are in states of plane stress 
It will be later explained under what conditions this is true. 
Then the stress-strain relations of a single lamina referred to 
its material coordinate system are written in the forms: 


F (k)_ q (k) 
e a8 a8y6 



(a) 


e (k)= S (k) 


c (k) 


(b) 


(3.1.2) 
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where (3.1.2a) is in tensor notation with range of subscripts 
1, 2 and (3.1.2b) is in matrix notation. It should be noted 
that (3.1.2) represent the stress-strain relations (2.1.9 - 
. 10) , i . e . , 


. 00 . 
'll ' 


y OO 

11 

,(k) 


(k) 


, (k) 


,(k) 

22 


. (k) 

'22 " 


. (k) 
12 


,00 

A 


TOO 

J A 

,(k) 

J 12 


2G 


(k) 

A 


,0O 

11 


(k) 

22 


Ok) 


(3.1.3) 


Let a laminate of rectangular 
by a uniform edge stress: 


11 (+a. 

x 2 ) 

12 ( ± a ' 

*2> * °12 

12 (X 1' 

- bl ■ °12- 

22 (x l' 

±b) - a 2 J 


form. Figure 8, be loaded 


(3.1.4) 


The elasticity solution of the laminate must satisfy the 
following requirements: 

(a) Equilibrium of stresses, 

(b) Traction continuity at laminae interfaces, 

(c) Boundary conditions (3.1.4), and 

(d) Displacement continuity at laminae interfaces. 

It is assumed that the stresses in any lamina are con- 
stant, but different in the different laminae. The condition 
(a) is satisfied within any lamina. Since the assumed lamina 
stresses are plane there are no traction components on laminae 
interfaces. Therefore (b) is satisfied. 

The boundary conditions (3.1.4) cannot be strictly satis- 
fied in each lamina but only in an average sense. To do this 
(k) 

lamina stresses o' referred to lamina material coordinates 


22 . 


are transformed to laminate axes. The stresses in the kth 

(k) 

lamina referred to laminate axes are denoted a a g‘ The 
transformation is given by: 

(k) (k) 2. 

a ll' a ll cos e k 

=„(k) _2_ 


(k) °l2- (a ll la i^ sin9 k 

or in matrix notation: 
(k) c = fl (k) 0 (k) 


f (k) 

22 

Sin \ 

- »£’ 

cos0. 

k 

sin 

6 k • 

(k) 

22 

2 

cos 6 k 

+ ^ 

COS0, 

k 

sin 

0 (3.1.5) 

°k 

cos9, 

k 

+ a (k) 
12 

o 

(cos 0, 
k 

- sin 2 

V 



(3.1.6) 


Let the edges of the laminate be loaded by constant forces 
per unit length T^, T 2 2 ' T i2 and def ^ ne the stresses (3.1.4) 
as edge averages over the laminate thickness h: 


'll 


22 


= T ll /h 
* T 22 /h 


12 


“ T 12 /h 


(3.1.7) 


Equilibrium requires that: 


K 

Z 

k=l 



o 


11 


K 

Z 

k=l 



o 


22 


K 

Z 

k=l 



12 


(3.1.8) 


where K is the number of laminae. Written in terms of stresses 
(k) 

o R using (3.1.6), we have: 


£ e (k) o (k) - a° 
k=l 


(3.1.9) 


where o° denotes the stresses o° at the edges. 
— ap 
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Replacement of the boundary conditions (3.1.4) by (3.1.6) 

is an approximation of Saint Venant type. Thus, there must be 

expected edge perturbations (among them interlaminar shear) on 

the stresses predicted by laminate theory. 

Equations (3.1.8) are three equations for the 3K stresses 

o^, in the laminae. There are needed an additional 

aB aB aB 

3(K-1) equations which are provided by displacement continuity 
at lamina interfaces, requirement (d) . 

Since the stresses in each laminae are by hypothesis uni- 
form, so are the strains. Therefore, displacement continuity 
is ensured if the lamine strains in adjacent laminae, referred 
to laminae coordinate system are the same. Thus: 




(k+1) 

£ 11 

(k+1) 

22 



(k+1) 

12 


k=l ,2 .k 


(3.1.10) 


Equations (3.1.10) are the additional required 3(K-1) equations. 

(k) 

They will be written in terms of laminae stresses o ' referred 
- ' aB 

to laminae material axes. To do this it is noted that: 


= £< k > £ < k > 

which is just a transformation of (3.1.6). From (3.1.2b): 


(k) e = e (k) s (k) a (k) 


(3.1.11) 


and inserting the last result in (3.1.10): 

e (k) s (k) a (k) = 9 <k+l) s (k+l) a (k+l) k=l , 2 k (3>1 . 12) 


Equations (3.1.9) and (3.1.12) are 3K linear equations for the 
3K stresses in an elastic laminate, with K layers. 

It should be carefully noted that the analysis given above 
is based on plane stress conditions in individual laminae. This 
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is a valid assumption if: 

(a) The loads on the laminate are statically equiva- 
lent to in-plane forces (membrane forces) and pro- 
duce neither bending nor twisting moments , and 

(b) The laminate has a certain stacking sequence of 
laminae which defines a so called balanced or 
symmetric laminate . 

This stacking sequence is an arrangement in which the 
laminate has a middle plane of geometrical and of material 
symmetry. The laminae are arranged in paris with respect to 
the plane of symmetry. The laminae of such pair have equal 
thicknesses, same distances from middle plane, and are of 
the same material with same angles of reinforcement. 

In a non-symmetric laminate application of membrane 
forces will in general produce bending and twisting of lam- 
inae and thus a plane state of stress will not be realized. 

The symmetric laminate is, however, sufficiently versatile 
to cover most cases of practical interest. 

Let it now be assumed that the laminate is inelastic but 
still fulfills the conditions of symmetry and pure membrane 
loading. In this case the only equations which necessarily 
change in the preceding development are the stress-strain re- 
lations of the laminae, (3.1.2), which must be replaced by in- 
elastic laminae stress-strain relations are given by (2.1.7) 
where the compliances are now functions of the stresses. These 
compliances now replace the elastic compliances in (3.1.2) 
which thus become non-linear. 

It is convenient for later purposes to rewrite (3.1.2) in 
the inelastic case in different form. To do this the strains 

are first split into elastic strains (2.1.9) and inelas- 
ats 

tic strains (2.1.11). Preceeding to (3.1.12) this equation 
assumes the form: 

0 (k+l) s l(k+l) a (k+l) _ 0 (k) s l(k) o (k) = 

_ fl (k+l) s ll(k+l) 0 (k+l) + e (k) s ll(k) 0 (k) (3.!.! 


k = 1, 2 


k 
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where 


gl (k) - elastic compliance matrix of kth layer 

gll(k)- inelastic part of compliance matrix of kth layer 

s (k) = s l(k) + s ll(k) (3.1.14) 

Equations (3.1.13) are now written out in component form 
with notation (2.1.10), (2.1.12) for compliances: 
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in ^ k U^ ^ k *Rns8. sin8. (3.1.15) 
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(k+1) ,„1 (k+1) . 2 . , _1 (k+1) 2. 

°11 (S 11 sin 9 k+l S 12 cos 9 k+l 

^ (k+1) , c i (k+1) . 2, ^ „1 (k+1) 2 a , 

+a 22 (S 12 Sln 9 k+l + S 22 COS e k+l> 
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12 12 44 k k 
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(3.1.16) 
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(3.1.17) 
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To these must be adjoined equations (3.1.9) which are written 
here in components : 

E (a|^os 2 0 k +O 22 ^ sin 2 0 k ~ 2 a| 2 sinO^) t fc = h (a) 

k=l 

k 

E (o j^sin 2 0 k +a22^cos 2 0 k +2a^2^ cos0 k sin0 k ) t fc = a° 2 h (b) 
k 

£ (a^ ~ a 22 ) cos0 k sin0 k +a l2^ cos2e k -s;i ' n2e k* fc k = °12 h 

(3.1.18) 

We now consider special cases of interest. In the first 
case the inelastic laminae strains have the form (2.1.13). 

Then the right side of (3.1.15-.17) simplifies by setting: 
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b 44 “ b 44 °22 ' °12 


(3.1.19) 


Once the stresses in the laminae have been obtained the 
strains in the laminae, referred to laminate axes, are deter- 
mined from (3.1.11). Since the strains in all laminae are 
the same when referred to the laminate coordinate system, these 
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are also the average laminate strains and thus determine the 
inealstic response of the laminate. 

In the simplest case the lamina material is assumed to 
be inelastic in shear only. In that event we have in addition 
to (3.1.19) : 




Oct!) = Q 


(3.1.20) 


and for Ramberg-Osgood presentation of inelastic part of shear 
compliance: 
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(3.1.21) 


In Ramberg-Osgood representation (2.2.1) the inelastic parts 
of the compliances assume forms such as : 
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3. 2 Method of Solution 

The equations which define the laminae stresses are (3.1. 
and (3.1.13) in condensed form, or equivalently, (3.1.15- 3.1 
(3.1.18) in full form. To explain the solution method it is 
simpler to write in terms of the condensed form. 

Define the matrices : 

L l(k+l) =e (k+l) s l(k+l) 

Lll(k + 1)- e ( k + l) sll (k+ l) ' (3.2 

L l(k) = 0 (k) gl(k) 

L ll(k) = Q (k) gll(k) 

Then equs. (3.1.13) assume the form: 

1,1 (k + D g (k+1) _^1 (k)^(k) = _ L ll(k+l) a (k+l) +L ll(k) 0 (k) (3.2 


to which are adjoined equs. (3.1.9) which are here rewritten: 


I e (k) o (k) = a 0 (3.2 

k=l 


The equations may be solved numerically by an iteration 

method which proceeds as follows: Consider equs. (3. 2. 2-3) 

with the right side of (3.2.2) zero. This defines a set of 
(k) 

stresses o„ given by: 


l(k+l) (k+1) T 1 (k) (k) 

_ 2.0 - t 2.0 = 


(a) 


k = 1,2 


. k=l 


(3.2. 


K 

Z 

k=l 


0 (k) a „ (k > = a < 


Since (3.2.4a) contains only elastic compliances S 


(b) 

, (k) 


it is 


seen that the equations are linear and define the stresses in 


an elastic laminate. Now insert the stresses a, 
right side of (3.2.2) and define the stresses 


(k) 

(k) 


into the 
by: 


9) 

17) , 


1) 


2 ) 


3) 
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l(k+l)„ (k+i) _ L l(k)^ (k) = 


2-1 


°1 


+L 11(k) [ Oo (k) ] a 0 (k) 


L H (k+1) [Q (k+I)j a (k+1) (a) 

(3.2.5) 


* e< k ’ o ,»> 

k=l 1 


o" 


(b) 


Equs. (3.2.5) defines (hopefully) a new approximation cr^ k ^ 
which is the solution of a set of linear equations. The stresses 
in square brackets in the right side of (3.2.5) are to emphasize 
the stress dependence of the non-linear parts of the compliances. 

The procedure just initiated can be repeated indefinitely. 

In general : 


L 1 (k + Dg(k+l)_^l (k)(j (k) = _l11(*+1) [o (k+Dj g (k+1) 

&+1 ~ £+1 “ ~i “fc 

+L 11 (k) [c (k) ] a (k) 

£ ~£ 


K 

L 

k=l 


e (k) a (k) 

£+1 


a° 


(3.2.6) 


This iteration procedure is quite easy to carry out with 
aid of a computer. It replaces the solution of a set of non- 
linear equations by solution of a sequence of linear equations, 
provided of course, that convergence is obtained. 

It should be noted that the first iteration step does not 
necessarily have to start with equs. (3.2.4a), i.e., with zero 

i k \ 

right side of (3.2.2). Any stresses o<> which fulfill (3.2.4b) 

can be used to start the iteration with (3.2.5) and continuing 

with the general iteration relation (3.2.6). 

It is desired to obtain a laminate solution for only one 

load system o° then it would seem most logical to start with 

(3.2.4). But suppose there is a sequence of loadings Ao ° , 

2Ao°...nAo . Suppose that a solution for (n-1) Ao° has been 
0 

obtained and that a solution for nAo° is desired. One possi- 
bility is to multiply all stresses due to the load (n-1) A o° 
by the factor n/(n-l). The stresses thus obtained certainly 
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also satisfy (3.2.6b) because of the linearity of these equa- 
tions. They will generally be reasonable starting values 
(k) 

a for the iteration. 

O 

This method of iteration to obtain a solution was found 
to work well for many sample problems; however, there were 
cases in which the solution did not converge. Attempts to modi- 
fy the recurrence relations to overcome this problem met with 
only partial success. Thus, ah alternate procedure for solu- 
tion was defined. The solution was obtained by application of 
the Newton-Raphson method. 

The set of 3K nonlinear equations represented by equs. 

(3. 2. 2-. 3) may be presented in the form: 

F n (o ij k) = 0 n = 1,2 . . . 3K (3.2.7) 


The function F^ is expanded in a Taylor series about an arbi- 
trary set of initial stresses which may be taken as the solu- 
tions of the elasticity problem. Considering only two terms 
of the series , it is found that 


F. 

1 


F. 

1 


+ 



3a 


mn 



mn 


0 


or 


a . . 

13 


O 



(3.2.8) 


(3.2.9) 


where a^j 
solution ' 


is the corrected solution obtained from the assumed 




Using aV^. as the initial guess, the process 


is repeated until the result is obtained within a desired ac- 
curacy. A recurrence form of equation (3.2.9) to obtain the 
stresses at t+1 cycle from t cycle can be constructed as 
follows : 
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(3.2.10) 



After the stresses a. . are obtained for all layers of the 
laminate, strains for any layer k in terms of laminae axes 
can be computed using equs . (3.1.3). Strains in terms of 

the laminate axes can be obtained using the strain trans- 
formation law. 


This analysis has been developed into an efficient 
computer program. A description of the program including 
a listing, is presented in Appendix E. 
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3.3 Numerical Results 


The computer program which has been developed under the 
present study has been utilized in the analysis of a variety 
of different composite laminates. The initial studies using 
the computerized analysis were directed at presenting a compari- 
son between the results of the present analysis and those of 
previous analyses, notably that of Ref. 9. (The present results 
were also compared to available experimental data, primarily 
those of Ref. 6 which had also been used for comparison with the 
analytical results in Ref. 9.) The objective of this phase 
of the numerical study was to determine whether the present 
results, which can be obtained with minimal computer usage, 
compare well with those of the more exact and complex analyti- 
cal results in Ref. 9. The results of this comparison are highly 
encouraging, as will be shown below, and support the utilization 
of the present analysis as an efficient design tool. 

In the second phase of the design numerical studies, con- 
sideration was given to examining the sensitivity of laminate 
results to individual properties of the layers. These para- 
metric studies are presented for several classes of typical 
laminates . 

A series of laminates of boron/epoxy composites for which 
experimental data had been obtained in Ref. 6 were examined 
analytically in Ref. 9. In Figures 9 to 15, results of the pre- 
sent analytical method are added to the comparison of experimen- 
tal results of [6] and analytical results of [9]. For example, 
in Fig. 9, the experimental stress-strain curve for a 0-90 boron/ 
epoxy laminate is compared to the analytical results obtained 
in Ref. 9 and in the present analysis. Both analytical re- 
sults coincide; both show slightly less inelastic strain than 
the experiment. The solid point on the curve indicates the 
stress level at which fiber fracture is computed to occur in 
one of the layers of the laminate. 

The shear stress-strain curve used in the present analysis 
was the best fit Ramberg-Osgood curve having an exponent n=3. 
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The values of modulus and yield stress obtained from the least 
squares fit are shown on the figure. A similar result is shown 
for the unidirectional tension +45° laminate in Fig. 10. Here 
it is seen that the two analytical curves are similar, although 
the agreement is not as close as in Fig. 9. Experimental data 
reflect a substantially higher degree of inelasticity than 
either analytical result. The present analysis shows a higher 
degree of inelastic strain at the higher stress level than that 
of Ref. 9. However, the reverse is true in the comparison of 
the two analytical results shown in Fig. 11 for a +30° laminate. 
The present results were obtained with a linear stress-strain 
curve in the transverse direction within each of the layers. 

The computations were made in this fashion because the transverse 
stress-strain curve of Ref. 9 does not show a significant degree 
of inelasticity. 

Figure 12 presents results for the case of a quasi-isotrop- 
ic laminate (0/+45/90) of boron/epoxy. Both the present result 
and that of Ref. 9 show a relatively insignificant amount of 
inelasticity. Again, the experimental data show a greater 
inelastic effect. Here the predicted failure strain level is 
in good agreement with the experimental failure strain level; 
however, there is a significant difference in the failure 
stress level. A similar result is presented in Fig. 13 for 
the quasi-isotropic laminate formed from the 0/+60 0 configura- 
tion. 

Computations performed for the present study for laminates 
having fibers in several directions, including the loading direc- 
tion, for a simple unidirectional load have shown a relatively 
small amount of inelastic strain. Another example of this is 
presented in Fig. 14 for a 0/+45 0 laminate. Here, however, the 
agreement of all the analytical methods and the experimental 
method is very good. 

The final comparison taken from Ref. 9 is presented in 
Fig. 15 for a laminate having fibers in three different direc- 
tions and a tensile load applied at some intermediate angle. 

The present analysis agrees reasonably well with the results 
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Ref. 9. The discrepancy between the failure load predicted 
on the basis of fiber failure and the experimentally ob- 
served failure stress is quite substantial. It is possible 
that fialure in laminate of this type caould result from 
shearing or transverse stresses within the individual layers, 
and thus, not be a result of tension in the fiber failure. 

This mode of failure has not been treated in the present 
computer program. The mode of failure observed experimentally 
is not known to the authors. 

The experimentally measured response of a multidirectional 
laminate to an applied shear stress has been reported in Ref. 
13. Comparison of the experimental result with the theory of 
Ref. 9 was presented in Ref. 14. Computations for this case, 
made using the present analysis and the prior analytical re- 
sult (Ref. 14), are compared to the experimental result in 
Fig. 16. Again, correlation between the two analytical re- 
sults is good, agreement between analytical and theoretical 
results is reasonably good with the experimental observation 
showing higher inelastic strains and lower tangent shear moduli 
at the very high stress levels. 

The conclusion of these comparisons with analytical and 
experimental data seem to justify the adoption of the present 
computer program as a useful engineering tool for the design 
and analysis of composite laminates. However, it appears that 
further study of the failure region is required. 

Parametric study of the influence of various laminate 
geometric and mechanical properties has also been explored. 

Fig. 17 shows the results obtained for a 0/+45 0 laminate in- 
dicating that the inelastic response in the transverse direc- 
tion can become significant at higher stress levels. Failure 
due to fiber fracture under a transverse stress applied to 
the laminate occurs at strain levels larger than those plotted 
in Fig. 17. In the quasi-isotropic laminate having four fiber 
directions, (0/+45/90) the degree of inelasticity in the longi- 
tudinal and transverse directions is of course the same and is 
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in both cases very small. It is to be anticipated, on an 
intuitive basis, that the maximum degree of inelastic response 
would be observed for a stress applied midway between two of 
the fiber directions on this quasi-isotropic laminate. The 
stress-strain curve for this latter case is also shown in Fig. 

18. Although the inelastic strains for this case are not 
significant there is a large difference in the predicted 
failure stress levels based on stress in the fiber direction 
for the two cases. It is worthwhile to emphasize that the 
quasi-isotropic laminate need not be isotropic in its strength 
characteristics . 

Because of the directional strength characteristics inter- 
esting effects may be expected for combined stress cases. Some 
results of the exploration of this question are presented in 
Fig. 19 where the four direction quasi-isotropic laminate is 
subjected to combined stress state with respect to a 22-1/2° 
axis of symmetry. This laminate shows high strength under both 
the unidirectional load and shear load by itself. The combined 
stress case for equal values of applied shear stress and axial 
stress results in fiber failure, and- therefore, laminate failure, 
at a substantially lower stress. The stress-strain curve prior 
to failure is not affected significantly by the presence of com- 
bined stress. The quasi-isotropic laminate having fibers in 
three directions (0/+60) is examined in Fig. 20. The sensitivity 
of this laminate to the Ramberg-Osgood parameters for the indi- 
vidual ply had little effect upon the stress-strain result. In- 
deed as an extreme example of this variation all lamiates stiff- 
nesses except the axial stiffness were equal to zero. Enforcement 
of the Kirchhof f-Love plate assumptions for this case results 
in the so-called netting analysis. The response for this net- 
ting case, which is linear, is shown by the dashed curve in 
Fig. 20. Even with this extreme assumption, matrix inelasticity 
does not introduce a significant amount of inelastic strain. 
Experimental data for comparison with this result are not easily 
available, however Ref. 17 does present a stress-strain curve 

for this case which shows a transverse failure stress for the 
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quasi-isotropic 0/+60® laminate which is about 60% of the 
failure stress in the axial direction. Also, the inelastic 
strain at failure is approximately 30% larger than the elas- 
tic strain associated with the failure stress level. The net- 
ting analysis result presented here suggests that in order to 
obtain such a strain, one might have to consider that the 
axial stiffness, either in tension, compression or both; or 
that other effects not considered in the conventional lamin- 
ate analysis, such as interlaminar or transverse shear de- 
formations, might contribute significantly to the overall 
laminate deformation. 

The influence of the characteristic stress levels for 
transverse stress and axial shear of the unidirectional layer 
of a boron/epoxy material is examined in Fig. 21. The measure 
of this effect is taken to be the influence upon the stress- 
strain curve for the unidirectional tension of +30° laminate. 
The strong sensitivity to the characteristic axial shear 
stress Ty and the relative insensitivity to the transverse 
characteristic stress a ^ for the R-0 representations is 
illustrated in the figure. . A similar comparison made for a 
boron/aluminum laminate of the same geometry subjected to 
uniaxial applied stress is shown in Fig. 22. Similar sensi- 
tivities are observed for this case. Boron/aluminum laminate 
response under transverse applied stress with the same values 
of the Ramberg-Osgood parameters is shown in Fig. 23. Here 
the fiber failure criterion did not come into play and thus 
the computations were extended to rather large strains in 
matrix. It is clear, that for this case, the failure criterion 
based on other stress-strain components is required. The exam- 
ination of the computer print-out permits one to terminate 
the stress-strain curves at some stress level prior to fiber 
fracture depending upon the choice of the failure criterion. 
This can be done rather readily. The choice of the failure 
criterion is discussed in Appendix D. 

The lamina properties for boron/aluminum are used to 
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analyze a 0°/+30° laminate under combined loading. These re- 
sults are shown in Fig. 24. Axial stress-strain curves are 
presented for varying ratios of axial shear stress to axial 
tensile stress. 
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CONCLUDING REMARKS 


Current approaches to the definition of design allowable 
stress for advanced fiber composite laminates are based upon 
the utilization of extremely conservative criteria. These limit 
the laminate to stress levels below which no significant damage 
of any kind occurs. The utilization of overly conservative de- 
sign criteria can negate much of the potential for effective de- 
sign utilizing advanced composite materials. The heterogeneous 
nature of these materials is such that a variety of possible 
damage modes exist. Thus, matrix cracking or yielding, fiber 
fracture, debonding, and other inelastic effects can all occur 
in local regions at relatively low average stress levels. These 
nonuniform and nonlinear effects greatly complicate the problem 
of establishing reliable design allowables. In the present pro- 
gram, the problem of nonlinear laminate behavior resulting from 
nonlinearities in the behavior of the matrix material was studied. 
The objective of the program was to develop an understanding of 
the inelastic behavior of composite laminates and to develop a 
computer program which will be used as an engineering tool in 
the design of fiber composite laminated structures. 

The method of approach utilized herein was to adopt a Ram- 
berg-Osgood representation of the nonlinear stress-strain be- 
havior and to utilize deformation theory as an adequate repre- 
sentation of the material nonlinearities. The problem was viewed 
on two levels. First, the relationship between the constituent 
properties and the stress-strain response of a unidirectional 
fiber composite material was studied. For this problem, the 
primary attention herein was directed toward the axial shear be- 
havior, in as much as experimental data had indicated that it is 
this type of load which results in the most significant nonlineari 
ties in material behavior. For this case, an expression was estab 
lished relating the composite average-stress/average-strain curve 
to the fiber moduli and the matrix nonlinear stress-strain curve. 
This expression, which was developed as a lower bound, was found 
to give good agreement with the more exact results obtained by 
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applying incremental plasticity theory and using a numerical 
finite element analysis to the assessment of the material be- 
havior (Ref. 9). 

The second level of approach treats the interelationship 
between the properties of the unidirectional layers and those 
of the laminate. For this case, one may consider that the 
starting point is a nonlinear stress-strain curve for trans- 
verse stress, and for axial shear stress, alone, and a linear 
stress-strain relation for stress in the fiber direction. The 
nonlinear lamina stress-strain curves can be modeled by proper 
selection of the Ramberg-Osgood parameters. 

In the present study, unlike other formulations an inter- 
action expression was formulated to account for simultaneous 
application of axial shear and transverse stress. A laminate 
having an arbitrary number of oriented layers, and subjected 
to a general state of membrane stress, was treated. The results 
of this analysis were programmed into an efficient computer 
routine for numerical evaluation of arbitrary laminates. Results 
obtained show good agreement with those of previous complex 
numerical methods utilizing incremental plasticity theory. 

Certain limitations connected with this program should 
also be discussed. First, deformation type stress-strain re- 
lations have been used; hence, it is implicit in this result 
that the stress and strain values obtained for any given set 
of loads are functions only of those loads and not of the 
loading history. On the other hand, if points are computed for 
intermediate values of loads, following different load paths, 
then different intermediate conditions will be obtained. Thus, 
the question is raised as to what is the accuracy of the results 
obtained for paths which do not yield proportional loading. It 
is known that for local proportional loading, the deformation 
theory result is the same as that for the incremental theory. 

In the laminate, local proportional loading does not exist, 
in general, even when the external loading is proportional. How- 
ever, the assumption is made that the deformation theory will 
yield an approximation which is satisfactory to generate a 
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rational engineering tool. This can only be assessed by com- 
parison with an exact analysis, or since this does not exist 
for the case of arbitrary loading paths, perhaps by comparison 
with experimental data. 

Comparisons of the present results with experimental data 
tend to show moderately good agreement. There are, however, 
cases in which experimental results show a higher degree of 
inelastic strain than predicted by the present analysis. These 
experimental data are quite limited and may be insufficient 
for drawing conclusions in this regard. 

The question of failure criteria incorporated into the pre- 
sent analysis required further consideration. The present 
analysis obtains more accurate representations of the stress 
components in the individual layers than have been obtained 
from elastic analyses. Hence, the use of these stress components 
in any failure criteria should represent an improvement in 
failure prediction 

In addition to a description of the methods of analysis, 
and of the numerical comparisons which have been carried out, 
the present report also presents a description of the computer 
program for study of nonlinear behavior of laminates in suffi- 
cient detail to permit the utilization of this program by 
others . 
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APPENDIX A 


SYMMETRY SIMPLIFICATION OF NON-LINEAR STRESS-STRAIN RELATIONS 

The most general inelastic stress-strain relations of the 
deformation type are of the form 

e ij = S ijkl °kl (1) 


where S. .. ,are functions of the stresses. Let it be assumed 
i]kl 

that the material is transversely isotropic with x^ axis of 

symmetry. Any rotation about x^^ changes . and ck . into e' . . 

and a ... Then the condition of transverse isotropy demands 
xj 

that 




= S. ., , o , . 
x]kl kl 


( 2 ) 


where in (1) and (2) are the same. To fulfill this last 

requirement it is necessary that be functions of stresses 

only through stress expressions which are invariant for rota- 
tions about the x^ axis. There are five such invariants and 
they are given by, [15] 


*1 = °11 


I 2 ~ °22 +0 33 1 3 °12 + °13 


(3) 


I 4 = 1 /2«, 22 -<,3 3 ) 2 + 2 023 2 


I 5 “ 1/2<0 22 _0 33 ) t0 12 °13 )+2 o 12°X3°23 
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Thus 


S ijkl S ijkl (I 1' Z 2' I 3' I 4' I 5 ) 


(4) 


It follows that for rotations around the x^ axis of 
symmetry the behave as constants. Consequently, the 

symmetry reduction of (1) to transverse isotropy is just 
as in elasticity. 

The reduction may be performed in following fashion: For 

rotation of angle 0 about the x, axis, the stress tensor a . 

i 1 13 

transforms into a ^ in the following fashion 


11 


°11 








22 

= 

1/2 

(0 22 

+ a 33 ) + 1/2 

(a 22 

C 33 ) cos 

20 

+ °2 3 s ^- n 

20 

33 

= 

1/2 

( °22 

+ °33> - 1/2 

(a 22 

a 33 ) cos 

29 

- a 23 sin 

29 

23 

= 

1/2 

( °33 

- a 22 )sin 29 

+ °23 

cos 29 




12 

= 

°12 

cos 

0 + a 13 sin 9 




(5) 


13 

= 

"°12 

sin 

0 + a^ 3 cos 0 







The same transformation relations obviously also hold for 
strains. If the transformed stresses and strains are intro- 
duced into (2) then coefficients of cos 29, sin 29, cos 0 and 
sin 0 and remaining terms independent of 9 must be equal. 

These equalities result in relations among the various compo- 
nents which reduce the stress-strain law to the form (2.1.4- 5) 
from Chapter 2 of this report. (Average stresses and strains 
appear in the latter but this obviously makes no differences 
in the derivation.) 
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APPENDIX B 


PLANE STRESS-STRAIN RELATIONS OF FIBER REINFORCED MATERIAL IN 
GENERALIZED RAMBERG-OSGOOD FORM 


The purpose of the present appendix is to arrive at 
equs. (2.2.7). For convenience in writing, overbars on 
stresses and strains will be omitted. 

The present development is guided by isotropic J 2 theory 
for deformation type plastic stress-strain relations. The basic 
assumption of this theory in the isotropic case is that the 
plastic strains have the form 


i: 


" f(J 2 )s ij 


(1) 


where 


s . . 


is the stress deviator and 


J 2 = 1/2 s ij s ij 


( 2 ) 


is its second invariant. 

It is instructive to work out the form of (1) for Ramberg- 
Osgood type stress-strain relations. Suppose that in pure shear 
the stress-strain relation is 


12 


12 


a n-1 
[1 + (r~~) 1 

2G y 

Now in pure shear it follows from (2) that 
_ 2 
J 2 12 

Therefore (3) can be written in the form 

/jT n-1 
[1 + (- — — ) 1 

y 


(3) 


12 


“12 

2G 


(4) 


which is in the form(l) . Consequently, in the general case of 
three dimensional stress and strain 


II 


e 


ij 



/JT n-1 

] 


y 


(5) 
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It should be emphasized that there is nothing fundamental 
about (1) . It is an assumption which states that the plastic 
strains can be represented by the stress deviator components 
multiplied by a function of a quadratic expression in the stress- 
es which is J 2 . ' The choice of J 2 for a quadratic expression 
is not arbitrary but may be arrived at by isotropy arguments. 

In an anisotropic material it may be assumed by general- 
ization that plastic strains are given by 

£ ij = S ij f (L > (6) 

Where L is some general quadratic function of the stresses. This 
assumption will form the basis of the present development. 
Consider the stress-strain relations (2.1.13). It 

M || 

is assumed that s__ and s„. functions of the most 

22 _ 44 

general quadratic form in a 22 and 
Thus 


II 

22 

II 

S 22 

(ao 22 

+ Bo 22 

°12 + C0 12 ] 


II 

44 

'll 

cn 

45 * = 
4 ^ 

(a? 22 2 

+ B ° 22 

°12 + Ca 12 * 

(7) 


It should be noted that the material reacts in same fashion 
to positive or negative shear stress, therefore also in same 
fashion to some cj 22 together with positive or negative shear 
stress. However, the middle term in the quadratic changes 
sign with shear stress. Therefore, this term should be omitted. 

Now rewrite (7) in form 


s 


22 • 




+ 



II 






( 8 ) 
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where f 22 and are nondimens ional functions and a and g 

have dimensions of reciprocal of stress. If o 12 =0 the first of 
(8) assumes the form 

f 22 (a c 22 ) (9) 


’22 


For one dimensional a 
strain relation (2.2.6a) 

a. 


22 


1_ . 22 . M-l 

E„ a ' 

T y 

which can be written as 


22 


22 


1 ° 22 2 

t O i 2 

T V 


from the Ramberg-Osgood stress- 


It follows from (8) and (10) that 

2 1 
a = — =■ 


(10) 


( 11 ) 


and the function of f 22 is determined as (M-l)/2 power. 

In similar fashion, when a^=0 , the second of (8) assumes 

the form 

s 


44 


2G T f 44 °12 2) 


( 12 ) 


From the Ramberg-Osgood relation (2.2.6b) for one dimen- 
sional a 12 

.. , o, 


44 


^ / 12. N-l 

2G_ 'x ’ 

T y 


which can be written as 
1 


44 


°1 2 2 Hll 

2GL- "r 1 * ] 2 

T y 


(13) 


It follows from (12) and (13) that 


3 2 = -4 


(14) 
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and the function f is determined as (N-l)/2 power. 
Consequently (8) now assumes the form 


s 


22 



°22 ^ 
Ko 21 * + 
y 



M-l 

2 


" = _J_ 

S 44 2G t 

Then (2.2.7) 


°22 ^ 
y 



follows from (15) 


M-l 

2 

and (2.1.13). 


( 15 ) 
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APPENDIX C 


1. EXTREMUM PRINCIPLES OF DEFORMATION THEORY OF PLASTICITY 
i. Principle of Minimum Potential Energy 


Let 


a ij C i jkl e kl 


( 1 . 1 ) 


where C. , are functions of the strains. The strain energy 
ijkl 

density is defined by the path dependent integral 
e 


W = / a .. (e)de . ■ 
t = 0 3 3 


( 1 . 2 ) 


where £ is a concise notation for E ij . The strain energy 
U £ of a body of volume V is defined by 


U e = / v W e dV 


( 1 . 3 ) 


Let the surface of the body be subjected to the boundary con- 
ditions 


u i (S) = u^ on S u 


( 1 . 4 ) 


T i (S) = T*. on S T 


and let the body forces vanish. The potential energy u 

I 

is defined by 


u p ■ / v w - dv - / s T 

Define an admissible displacement field {y (x) by 


u. = u° . on S 
i i u 


( 1 . 5 ) 


u i (x) continuous everywhere 


( 1 . 6 ) 


Associated with 5, are the strains ? derived from it by 

. 1 ID 

the usual relations. 
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£ 

Define W by 
e 

W £ = / a . . de . . (1.7) 

13 13 

e=0 

where 


5 ij = C i jkl } £ kl 


Define 


( 1 . 8 ) 


Up = / v W e dV - / s T°.u.dS 


(1.9) 


The principle of minimum potential energy for the present case 
then states that 


U > U 

p - p 


( 1 . 10 ) 


equality taking place if and only if 


u i = U i 


In the event that displacements are prescribed over the 
entire surface, the surface integral in (1.9) vanishes. Then 
the principle reduces to that of minimum strain energy 


U 


£ 



( 1 . 11 ) 


ii. Principle of Minimum Complementary Energy 
Let 


e . . 

1: 


“ S i jkl ( - ) °kl 


where S- . ate stress dependent compliances 
13 k! 


( 1 . 12 ) 
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Define the complementary energy density W° by the path 


dependent integral 

a 

W° = ! e . . do . . 

il il 

a=0 


(1.13) 


Let the surface of the body be subjected to the boundary con- 
ditions (14) and let the body forces vanish. The complementary 
energy is defined by 


U c = / v W u dv - /gJiU*. dS 


(1.14) 


Define an admissible stress field 5 . . by the following 
requirements 


a . . , .=0 
il 1 


T. = a . .n . 
i il 1 


continuous everywhere (1.15) 

T ± (S) = T°. on S T 

Define the complementary energy functional U c by 


0 - V* a dV - / T „• dS 

u 


where 


W ° = £ ij d5 ij 
5=0 


(1.16) 

(1.17) 


£ ij = S ijkl ( i } 5 kl 


Then the principle of minimum complementary energy states that 
U c >U C (1.18) 

equality occurring if and only if 


o . • = a . . 
il il 


If tractions are prescribed over the entire surface, S u =0 , 
the principle reduces to 

S' > u° 
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For proof of these principles see e.g. [16]. An interesting 
application to obtain approximate solutions has been given 
in 1 1 7 ] . 


iii. Specialization of the Principles to Axial Shear with 
Ramberg-Osgood Stress-Strain Relations 
In the case of axial shear of a uniaxially fiber reinforced 
material the only surviving stresses are 


0 12 T 2 


°13 = T 3 


( 1 . 20 ) 


where 1 indicates fiber direction. Denote the associated shear 
strains by 

( 1 . 21 ) 

£ 12 £ 2 £ 13 £ 3 

Then the generalized Ramberg-Osgood stress-strain relations. 
Appendix B, (5) assume in the present case the form 

t 2 t n_1 

e 2 - [ 1 + <h> ] 

y 


T 3 t n_ ^- 

’ 55 1 1 + ( h> 


( 1 . 22 ) 




T 2 + X 3 2 


In the present case 

a ij d£ ij = 2{T 2 d£ 2 +T 3 de 3 ) 


(1.23) 
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Inserting (1.22) into (1.23) and using the relation 


tdx = T 2 d x 2 +x 3 d t 3 


it is easily shown that 


a. . de 1 + n (— ) 

X] 13 G T y 


n -1 


] dx 


(1.24) 


To compute Iv as defined by (1.2) it is necessary to 
integrate (1.24) from zero to some state of strain e 2 , e^. 
But it should be noted that (1.24) is expressed in terms of 
the variable t only. Now t can be expressed in terms of 
strains in following fashion. Define 


t = / 2 2 

e 2 + 3 

It follows at once from (1.22) that 
x t n_1 

e = h [1 + 1 

This relation defines x as a function of 
quently, w e assumes the form 

e _ 1 T <e) n-1 

W G / x [1 + n (i-) ] dx 

o 

which is easily integrated to yield 

_ 2 0 n-1 

w e = ~ [l + 2 - n - (— ) 1 

2G 1 n+1 'x y ; J 


(1.25) 

(1.26) 
Conse- 


(1.27) 


x = x (e) 


0 

According to (1.3) the strain energy U is then given by 
the volume integral of (1.27). Note however that it is very 
difficult to express U £ in terms of strains since this requires 
the solution of (1.26) for x in terms of e. In general it is 
not possible to do this analytically. This places a severe 
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limitation on the use of the principle of minimum potential 
energy or of minimum strain-energy with Ramberg-Osgood stress- 
strain relations. 

Next we consider the principle of minimum complementary 
energy for axial shear. Since there are only shear stresses 
x 2 , and shear strains e^, the integrand in W ° , 

(1.13) , is given by 


£ ij da ij = 2(e 2 dT 2 +£ 3 dT 3 ) 


It follows from (1.22-.23) that (1.28) 


(1.28) 


is given by 


e. .da. . = £ [1 + (— ) 
13 13 G v T y ' 


n-1 


3 dt 


Integration of this expression from 0 to t 


yields 


,,0 1 r 

W = 2G [ 


1 + h£t <t-> 

y 


n-l 


3 


(1.29) 


Expression (1.29) now enters as the integral into the volume 
integral of U c , (1.14). 


We now examine the meaning of an admissible stress field 


2 ' 


in the present case. 


equation is 
3L 


3x, 


3x . 
3x- 


The only surviving equilibrium 

(1.30) 


= 0 


'2 3 

The traction components are 


T 1 " 


T 2 n 2 +T 3 n 3 


m ^ 
2 


T 2 n l 


T = 
3 


x 3 n l 


We shall be concerned with 
reinforced materials whose 
a surface n 1 =0. Therefore 
nent on such a surface is 


(1.31) 


cylindrical boundaries in fiber 
generator is in direction. On such 

the only surviving traction corapo- 



T 2 n 2 +T 3 n 3 


(1.32) 
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Consequently an admissible stress system t 3> must satisfy 

(1.30) and the value t° of t wherever prescribed on the 

n n 

boundary. 

The complementary energy functional (1.16) assumes the form 


W° = 


/ v W°dV 

- / T u° dS 
S 1 L 

(a) 

~ 2 

, u - n-1 


35 [1 

+ nil ^ 5 1 

y 

(b) 

/- 2 

~ 2 

(c) 

T 2 + 

T 3 


(1.33) 
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2. LOWER BOUND FOR AXIAL SHEAR MODULUS 

Consider a uniaxially reinforced lamina which is subjected 
to axial shear t c in the 1-2 plane on its boundary, fig. 5. 

By the average stress theorem, of Ref. 5. 

0 12 = To (2a) 

and all other average stresses vanish. 

By the average theorem of virtual work, of Ref. 5, 


/ e..do.. = e..d a.. 
y id id id ID 


( 2 . 2 ) 


Since the only nonvanishing average stress in the present 
case is (2.1) we have 

t. .do. . = 2~e dx 0 (2.3) 

ID ID 12 

The complementary energy of the body is given by (14) of 
Appendix A. The surface integral vanishes however in the pre- 
sent case since no displacements are prescribed on the 
boundary. Now 


u =; w°dv = / / 

L V 


0 _ 

e . . do . . dV 
T-D ID 

o=0 


(2.4) 


= 7 /e . .do . . dV 
ID ID 

o=0 


To _ 

f e 12 dT, 


The last equality following from (2.2, 3). 

s • 

By definition the effective secant modulus G is given by 


Hence 


(2.4) 
U C = V/ 


l Q 

G A<°12> _2G 1 

assumes the form 

T 0 T 0 dT 0 

° G l (T 0 ) 


(2.5) 


( 2 . 6 ) 


In order to find a bound on it will be necessary to 
find a bound on (2.6) by use of the principle of minimum com- 
plementary energy. 
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It is assumed that the fibers are infinitely rigid in com- 
parison to the matrix. Therefore at fiber/matrix interface 

n 1 = 0 ( 2 . 7 ) 

and the only contribution to the complementary energy is from 
the matrix. Thus, the surface integral in (1.33a) vanishes and 
it can be written as 

U = / W° dV (2.8) 

c V 

m 

where V is the matrix volume, 
m 

Furthermore, by (2.3.3) the actual stresses are functions 
of x 2 , x^ only. It is therefore natural to also choose admissible 
stresses as functions of x 2 , x^ . Thus W° in (1.33) becomes 
a function of x 2 , x^ only and therefore without loss of generality 
(1.33a) can be taken over unit length in fiber direction. Thus 
it can be written 

U c = f A w° (x 2 , x 3 ) dx 2 dx 3 (2.9) 


In order to construct an admissible stress system it is 
necessary to devise a geometrical model for a uniaxially rein- 
forced material . In past analyses of FRM two kinds of models 
have been successfully treated: Periodic arrays of identical 

circular fibers have been analyzed numerically with the aid of 
computers and the composite cylinder assemblage model has been 
treated analytically [1,5] yielding simple closed results. Since 
the present treatment is to be analytical the composite cylinder 
assemblage model will be used. A detailed description of the 
model has been given in [5]. Suffice it to say here that the model 
represents a cylindrical specimen of a fiber reinforced material 
as an assemblage of composite cylinders of different sizes which 
fill the space in the limit. In each composite cylinder the 
inner cylinder is a fiber and the outer shell is matrix material. 
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In all cylinders the ratios of fiber to matrix shell radius are 
the same, (figure 26). 

It is recalled that an admissible stress system must satisfy 
equilibrium and boundary conditions. An obvious possibility 
for such an admissible field are the stresses of the elastic 
solution since they certainly satisfy the required conditions. 
These stresses are the same in any composite cylinder of the 
assemblage and are given in cylindrical coordinates by (see [5]) 


rz 


= T r = t+c (1 + cos e 
r 


( 2 . 10 ) 


°ez T e “ i+c d sin 6 


0Z 0 l + c 1 ' 

where c is the volume fraction of fibers, a is the radius of any 
fiber and r,0are polar coordinates, fig. 26. 

Since t as expressed by (1.33c) is an invariant with 
respect to rotations about x, = z we have also 


2 2 ^2 

T = T + T. 

r 0 

Substituting (2.10) into (2.11) yields 


( 2 . 11 ) 


2 

T 


2 


= P 



cos 0) 


where 


( 2 . 12 ) 


P = 




(2.13) 


To simplify the analysis the exponent n in (1.22) will 
be assigned the value 


n = 3 


(2.14) 


It has been found that with this value of n, experimentally 
obtained shear stress-strain relations of epoxy can be quite 
accurately represented with proper choice of x . Recalling 
(1.33), (2.9) then assumes the form 

V C = 1G f A V t 1 + | <— (2 ‘ 15) 

m z y 
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where G is the matrix elastic shear modulus. Let the assemblage 
consist of K composite cylinder. Define U c for the kth comp- 
osite cylinder by 


U ^ = — / T 2 

C 2G A . 

mk 


[1 + i (J-) ] dA 
y 


(2.16) 


where A^ is the matrix area a^£r£bj, in the kth composite 
cylinder. Then 


U C = 


K ~ k 

z u c 


(2.17) 


k=l 


Since t 2 has been expressed in polar coordinates, (2.12), 
it is convenient to also evaluate (2.16) in the same coordinates. 
Using the variable p we have 



2 7T -2 

/ f 2 [1 + (£-) ] 
o y 


pdpde 


(2.18) 


where 

3 = V a k 


(2.19) 


which by construction has the same value in all composite cy- 
linders. Note also that the volume fraction of fibers c is 
given by 

, k. 2 _ 1 _ 

c(.) ^ (2.20) 

D k B 


Substituting 
tion we have 



(2.12) into 


(— 

l l+c 


(2.18) and carrying out the integra- 


, T 0 * 2 3+l0c' - 12c 2 -c 4 

( T y > TTT+cT 


(2.21) 
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is the area 


where (2.20) has been used. It is seen that TTb 2 ^ 
of the cross section of the kth composite cylinder and the 
parenthesis has the same value for all composite cylinders. 
Therefore, if (2.21) is inserted into (2.17) we find 


U C 2G T ° 2 [ 


1+ c 


+ (^) 


x 2 3+10 cr 12 i 


6 (1+ c) 


] 


( 2 . 22 ) 


Let (2.22) be written 


To ! du 

°C = A ' . A dx7 d 


(2.23) 


Without loss of generality (2.6) can be evaluated for unit 
height of cylindrical specimen. Thus 

U = A /° T ° dT ° (2.24) 

c ” '"•> 

Now introduce (2.23) and (2.24) into the minimum complementary 
inequality (1.18). Thus 


t 0 dU_ 

/ ri — - 

O l A dx 0 


To 

Gj(x 0 ) 


dx „ ^ 0 


Since the integral is positive for all values of 
the integrand must also be positive for all values of 
It follows that 


(2.25) 


To 
To • 


g a "•» i 


_AX a 

dU c /dx c 


J A(~) 


(2.26) 


where the extreme right denotes lower bound on the secant 
modulus g s . Substituting (2.22) into (2.26) and rearrang- 
ing we find the lower bound (2.3.9) of Chapter 2. 

There naturally arises the question of the establishment 
of an upper bound. The difficulties involved have been dis- 
cussed above: It is not in general possible to solve Ramberg- 

Osgood relations for stresses in terms of strains. It is 
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therefore not possible to analytically express the potential 
energy functional in terms of admissible strains. 

A possibility .to resolve the difficulty is to write in- 
elastic stress-strain relations of type (1.22) in the form 


2Ge 2 

[1 - 

<f-> 

a-!] 



y 




£ 

a-1 

3Ge 3 

(1 - 


] 


Y 



(2.27) 


where a and e are to be determined by curve fitting. The 
minus sign in the parenthesis is due to the fact that the 
stress-strain curve is below a straight line with the initial 
slope. 

It should be noted that (2.27) are not an inversion of 
(1.22). They are merely another form of approximation of 
actual stress-strain curves. 

In principle the representation (2.27) can now be used in 
conjunction with the principle of minimum potential energy 
to establish an upper bound on in same fashion as a lower 
bound has been established. It has however been found that in 
attempting to fit (2.27) to actual epoxy stress-strain curves 
a fractional exponent a was needed. This led to integrals of 
formidable difficulty in the evaluation of potential energy 
functionals. Therefore this approach has not been continued 
here. 
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APPENDIX D 


FAILURE OF NON-LINEAR LAMINATES 


It is expedient to separate the problem of the establish- 
ment of failure criteria of laminates into two separate problems: 

(a) Establishment of failure criteria for uniaxially 
fiber reinforced material, i.e., laminae. 

(b) Establishment of failure criteria of the laminate 
on the basis of laminae failure criteria. 

A great deal of wrok has been done on problem (a) . The 
problem has been approached in micro as well as macro-fashion. 

In micro-approach, it is attempted to predict failure on the 
basis of local analysis of the interior of the composite. Such 
an approach evidently encounters extreme difficulties. Although 
important work of fundamental nature has been done in this area, 
we shall not be concerned with it here since the work has not 
advanced to the stage of prediction of failure criteria under 
states of combined stress. 

In the macro-approach, a failure criterion is heuristically 
postulated as some function of pertinent state variables (gener- 
ally average stresses) which also contains undetermined para- 
meters. These parameters are then to be determined in terms of 
experimentally accessible information. 

We shall in the present discussion limit ourselves to states 
of plane stress. The simplest failure criterion is the so-called 
maximum stress criterion which states that failure occurs when 
either one of: stress in fiber direction, stress transverse to 
fibers, shear stress, reaches its critical value, these cri- 
tical values being the same whether or not the stresses act 
simultaneously. In symbols the criterion is: 

°11 = °A 


or 


22 


= a„ 


( 1 ) 


12 


T 


AT 
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or 


where 1 is fiber direction and 2 is the transverse direction. 

Generally, failure stresses o and o T are different in 
tension and compression. This is known as Bauschinger effect. 
There is evidently no Bauschinger effect for the shear stress. 
The simplest generalization of (1) to account for Bauschinger 
effect would be to assume as failure criterion: 


a 

11 

II 

Q 

> 

+ 

if 

a 

11 

°u 

°A 

if 

°11 

°22 

= 0 + 
T 

if 

°22 

°22 

°22 

if 

°22 

°12 

t at 

all 

°12 


( 2 ) 


whichever occurs first, where (+) and (-) superscripts denote 
failure stresses in tension and compression respectively. The 
main drawback of these simple criteria is in that they take 
no account of interaction effects. 

The most commonly used criterion which takes into account 
interaction is of quadratic form. For plane stress it has the 
form 


A ll°ll + A 22 °22 + A 12 a ll a 22 + A 44°12 = 1 (3) 

Here, products of shear stress with normal stress have been 
omitted since the material cannot distinguish between positive 
and negative shear stress. Therefore, odd powers (one, in this 
case) of shear stress cannot appear. 

Applying (3) to failure for stress in fiber direction 
alone, stress transverse to fiber direction alone, shear stress 
alone, in turn, it is seen at once that 
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( 4 ) 



The coefficient A^ 2 is troublesome since its determina- 
tion requires a failure experiment under combined stress. 
Several authors have proposed to use failure experiments on 
off-axis specimens under uniaxial stress for the determination 
of See e.g. [18] for discussion. 

The situation becomes more complicated if it is required 
to take into account Bauschinger effect, that is difference of 
failure stresses in tension and compression. One possibility 
to account for this effect is to assume that A^, A 22 assume 
different values for tension and compression. The situation 
regarding A 12 , however, becomes very awkward as it would have 
to assume four different values to account ,for four different 
possibilities of sign combination in biaxial stressing 
and 

It is also possible to add linear terms to (3) in which 
case it would assume the form: 

A ll a ll +A 22 a 22 + A 12 a ll°22 + A 44°12 + 


B l°ll + B 2°22 


1 


Such a device was suggested by Hoffman [19]. In this case it 
is possible to determine values of A^ , , A 22 , B 2 to account 

for different tensile and compressive uniaxial failure stresses 
in fiber direction and transverse to it. But the difficulty of 
assigning four different values to A^ 2 remains, unfortunately. 

In summary, the status of quadratic failure criteria has 
to date not been finalized. However, special versions of such 
criteria have been successfully fitted to experimental data. 

It is of importance to realize that in the fiber rein- 
forced materials used in practice failure predictions on the 
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basis of maximum stress criterion or quadratic failure cri- 
terion are not very different. This is due to the large 
ratios between strength in fiber direction and transverse 
and shear strengths and is easiest realized by considering 
the failure criteria as surfaces in cr.^, °22' °12 stress 
space. The maximum stress criterion is a very elongated 
rectangular parallelopiped while the quadratic failure cri- 
terion is an ellipsoid. For A ^ 2 =0 > Fig- 25 shows this 
schematically on a cut in the 0 -^, a 22 P^- ane - Thus it is 
seen that stress points on the two failure surfaces are close 
together for most parts of the surfaces. 

The situation would be entirely different for a material 
in which a , and a were of comparable magnitudes. 

A 

We shall now consider problem (b) i.e., the establishment 
of laminate failure criteria in terms of laminae failure cri- 
teria. The most conservative laminate failure criterion is to 
assume that once any lamina has failed the laminate has reached 
its ultimate load. There are cases of laminates in which all 
laminae would fail simultaneously and then this criterion would 
be justified. For example: a +0 laminate in which the exter- 
nal load direction bisects the angle between the fibers. 

In most cases, however, a certain group of laminae will 
fail first and failure of remaining groups would require fur- 
ther increase of load. Therefore a more realistic alternative 
is to determine the load at which the first laminae group fails. 
At this state, the further carrying capacity of the laminate 
may be assumed to be given by the remaining undamaged laminae. 
The increase in load which fails another group of laminae is 
then determined. This process is continued until failure of 
all laminae has taken place. 

Still another possibility is to assume that when a lamina 
has failed, certain of its stiffnesses reduce to zero. For ex- 
ample: suppose that a lamina or group of laminae has failed 
in shear. Such a failure implies a crack through the lamina in 
fiber direction. In that event, it is reasonable to assume 
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that the shear and transverse stiffnesses of the lamina are 
zero, but it still retains its stiffness in fiber direction. 
If, however, a lamina fails because of the stress in fiber 
direction the damage is so widespread that all of its stiff- 
nesses will be negligible. According to the type of failures 
encountered analysis is continued for the damaged laminate 
with the new stiffness rearrangement. This process is con- 
tinued until failure of all laminae has taken place. This 
method of analysis seems to be the most realistic but is 
also the most complicated. 

In almost all of the practical strength analyses of 
laminates in the literature, according to any of the methods 
outlined above, the stresses used for failure criteria have 
been determined on the basis of elastic laminate analysis. 
With the present inelastic laminate analysis, more realistic 
stresses are available in a better assessment of laminate 
failure loads. 
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APPENDIX E 


MSCtNOLIN COMPUTER PROGRAM 


1. General Description of the Program 

This is a computer program developed for the inelastic 
analysis of a laminate subject to any constant, arbitrary 
combination of in-plane loading. Details of the method of 
analysis and of the numerical solution, using the Newton- 
Raphson method, have been described in the body of this 
report. The essential features of the program are summari- 
zed below. 

The primary capability of MSC-NOLIN is to compute lam- 
inae properties when the laminate loads are defined. There 
is also a limited capability to work with constituent proper- 
ties, rather than laminae properties, as the input. Details of 
the input options are discussed subsequently. Basically, the 
inputs required are the stress-strain characteristics of the 
individual laminae for each of the three in-plane stress 
components applied separately. The stress-strain curves for 
transverse stress and for stress and for axial shear stress 
are defined by Ramberg-Osgood stress-strain curves. The 
parameters for these curves along with the laminae elastic 
constants are the required material property inputs. 

It has been observed that axial shear stresses in indivi- 
dual laminae are a major, perhaps the major, source of non- 
linearities in laminate response. Therefore, several additional 
options have been included in the MSC-NOLIN to accomodate more 
detailed characterization of shear response. First, the lam- 
inae shear stress-strain response may be input in tabular 
form and a least squares fit to the data ia automatically ob- 
tained for the R-0 yield stress (limited to the use of an 
exponent, n=3). Secondly, the matrix shear stress-strain curve 
can be input along with fiber elastic properties and the laminae 
shear stress-strain curve will be computed. In this latter case, 
the laminae elastic constants are also computed. 
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The input specifies one of two options for the determina- 
tion of the initial set of stresses to be used in the itera- 
tion at each value of applied load on the laminate. In one 
case the stresses found at one load are increased to the load 
for which the stresses were evaluated. In the other, and 
generally used option, the increment between the initial stresses 
used at the nth laminate load value and the actual stresses 
found for the (n-l)st load value bears the same relation to the 
ratio of those two load values as the similar relation computed 
at the previous load cycle, that is, 



V F 


(n-1) 
- a . . 

hi 

n-1 


(n-1) 




o 


(n-2) 

hi 


-;/ F 


' n-1 


n-2 


The program contains a number of controls to define: the 
size and number of steps of loading at which computations are 
made; the maximum number of iterations to be permitted in the 
numerical solution; the desired accuracy to be obtained in 
convergence; the criteria for divergence of the solution in 
the iterative process to avoid the use of unnecessary execu- 
tion time in the case of breakdown of the solution procedure. 

The program defines the failure of the laminate in a limited 
fashion, either on the basis of the maximum allowable stress 
in the fiber in tension or compression, or on the basis that 
the tangent modulus of the stress-strain curve of the laminate 
becomes less than a specified value. Failure due to shear or 
transverse stress are not included at this stage in the develop- 
ment of the program. 

2 . Input 

The main features of input in this program are the follow- 
ing : ; 

(a) Specify the number of laminates or problems to be 
solved; 
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(b) Define the geometrical properties of each layer; 

(c) Define either the material properties of each layer 
or the properties of its constituents; 

(d) Define either of the following for each layer: 

(i) yield stress in transverse direction and yield 
stress in shear; 

(ii) yield stress in transverse direction and a table 
of values defining shear stress-strain curve for 
the matrix plus a set of values of stresses to 
be used for the computation of yield stress in 
shear ; 

(e) Specify the type of Ramberg-Osgood relation to be used; 

(f) Define the loadings; and 

(g) Define the control parameters. 

A guide to the preparation of input data for this program 
is given in section 4 below. 


3 . Details of Output 

The output Coin be divided basically into two steps: 

(a) Output of Input Data: 

The first section of the output deals with the output 
of the input data. If the input is in the form of 
properties of constituents of the layer, it gives 
an output of the properties of the constituents first 
and then the computed value of the properties of the 
layer; otherwise, it gives output directly the prop- 
erties of the layer. 

(b) Output of Stresses and Strains: 

For each set of loading, the computer prints the 
following: 

(1) value of the load applied; 

(2) number of iterations . for convergence; 

(3) stresses for individual laminae with respect 

to principal elastic axes of the laminae; and 
strains for individual laminae in terms of both 
laminae and laminate axes. 69. 
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4. Input Details for MSC-NOLIN 
(1) Read (15) NSETS 


NSETS: number of problems 


(2) 

Read 

(15) 

LAY 


LAY: 

number of layers in this laminate analysis 

(3) 

Read 

(15) 

INP 


INP : 

Option for reading in material properties 


INP = 

1; 

read in material properties of individual 
laminae ; 


INP = 

2; 

compute properties of laminae from the 
properties of constituents. 

(4) 

(a) 

If INP = 1 


(i) 

Read 

(5D15.5) E n ,E 22 ' y 12 ,y 21 

(ii) 

Read 

(5D15.5) G 12 , SY, TY 

(iii) 

Read 

(D15 . 5,15) T, IANG 

(b) If 

INP = 

2 

(i) 

Read 

(4D15.5) EF, MUF, GF , VF 

(ii) 

Read 

(3D15.5) EM, MUM, GM 

(iii) 

Read 

(15) 12 


If I 

2=0; read in SY and TY 


(i) 

Read (2D15.5) SY, TY 


If I 

2 = 1; TY is to be computed 


(i) 

Read (5,1002) SYCE 


(ii) 

Read (215) NUMT 

NUMT = number of values in the table 


(iii) 

Read (5D15.5) TAU (J) , J=l, NUMT 
(Table of shear stress values of 
matrix read in) 


(iv) 

Read (5D15.5) GAM (J),J=1, NUMT 


(Table of shear strain values of 
matrix read in) 

(v) Read (5D15.5) SG12 (J) ; J=2,ll 

(Table of shear stress values of 
laminae read in) 
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(5) 


Read (5D15.5) XN, XM 


XN: exponent in nonlinear transverse stress- 

strain law; 

XM: exponent in nonlinear shear stress-strain 

law. 

(6) Read (5D15.5) SOU, S022, S012 

SOU: applied stress in X-direction 
S022: applied stress in Y-direction 
S012: shear stress in XY 

(7) Read (15, D15.5) KSGM, SMLT 

KSGM: total number of loading increments 

SMLT: ratio of load increment to the initial 

load. 

(8) Read (D15.5) STIFF 

STIFF: tangent modulus of stress-strain curve 

in terms of the laminate axes; specify 
a value of STIFF below which the program 
will not run. 

Read (ol5.5) SGR 

SGR: maximum allowable stress in the fiber in 

tension or compression 

(9) Read (I 5, 2D15.5) IT, EPS, UPBD 

IT: maximum number of iteration permitted in 

Newton-Raphson analysis 

EPS: convergence criteria; (ratio of values of 
two successive iterations should be less 
than EPS) 

UPBD: divergence criteria (solution will stop 
if ratio of two successive iterations is 
greater than ^+12^ 

(10) Read (I 5) INMT 

If INMT = 1, the program uses ratio of previous 
two solutions as the initial guess 
value iteration process; 

If INMT = 2, the program uses extrapolated value 

of previous two solutions proportioned 
on the basis of stress ratio as the 
initial guess. 
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Fiq. 2 - Nondimens ional Ramberg-Osgood stress-strain curves 
(e.a. . 2.2.1 with k =0.4). 
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^9* 3 - Transverse normal a-e curves for unidirectional Boron/ 
Epoxy (for normal loading and normal + shear loading) . 



Fig. 4 - Shear o-e curves for unidirectional Boron/Epoxy (for 
shear loading and normal + shear loading) . 
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Generalized Ramberg-Osgood stress-strain curve in axial 
shear for composite and matrix (n=3) . 







STRAIN 


Fig. 9 - Tensile stress-strain curve with [0/90.] B/Ep. 


82 . 



STRESS 



STRESS MN/m 



Fig. 11 - Tensile stress-strain curve with [+30] B/Ep. 


84 . 




STRAIN 


Fig. 12 - 0° Tensile stress-strain curves for [0/+45/90] B/Ep. 
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Fig. 13 - 0° Tensile stress-strain curves for [0/+60] B/Ep. 
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STRAIN 

Fig. 14 - o° Tensile stress-strain curves for [ 0 ^/+ 45 ] B/Ep. 
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Fig. 15 - 0°Tensile stress-strain .behavior of [65 3 /20/-70] B/Ep 
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Fig. 16 - Comparison of present results with experimental data 
for Glass/Epoxy 90/+18 tubes in torsion, said with 
the results of Ref. “9. 
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Fig. 17 - Stress-Strain curves for Axial, transverse and shear 
loading for [0°/+45°] laminate. 
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Fig. 19 - Combined stress effects on quasi-isotroDic [+22 . 5/+67. 5] 
laminate. ~ — 
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Fj-9- 20 - cj f feet of laminae inelasticity on transverse stress- 
strain curve for 0/+60 quasi-isotropic laminate. 
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Fig. 21 - Influence of laminae inelasticity upon axial tensile 
stress-strain curve of +30° Boron/Epoxy laminate. 
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Influence of laminae inelasticity upon axial tensile 
stress-strain curve of +30° Boron/Aluminum laminate. 
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Fig. 23 - 0° Tensile stress-strain curves for [0^7+4 5] B/Ep 
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Fig. 24 - Stress-strain curves of Boron/Aluminum [0/+3Q1 
laminate under combined loading. ~ 
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Fig. 25 - Comparison of failure criteria. 
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